# A “Watershed” Moment for Collaboration: 
# The Impacts of Legislation and Network Governance on Environmental Outcomes
# Yixin Liu, Jiho Kim, Tina Nabatchi 2025
# Public Administration Review
# Replication file for longitudinal regression analysis

library(dplyr)
library(spdep)
library(plm)
library(texreg)
library(purrr)
library(xtable)

get_summary_stats <- function(x) {
  tibble(
    N = length(x[!is.na(x)]),  # count non-NA values
    mean = mean(x, na.rm = TRUE),
    sd = sd(x, na.rm = TRUE),
    min = min(x, na.rm = TRUE),
    max = max(x, na.rm = TRUE)
  )
}

####Data####
# Setup working directory
WD = getwd()
setwd(WD)
df <- read.csv("df_for_longitudinal.csv")

####Table 2####
# Model 5
m5<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen) +
          lag(local) +
          lag(state) +
          lag(goals_ratio) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi1,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m5)

# Model 6
m6<-plm(wqi ~ 
          lag(cgr)*lag(goals_ratio) +
          lag(citizen) +
          lag(local) +
          lag(state) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi1,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m6)

# Model 7
m7<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen)*lag(goals_ratio) +
          lag(local) +
          lag(state) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi1,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m7)

# Model 8
m8<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen) +
          lag(local)*lag(goals_ratio) +
          lag(state) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi1,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m8)

# Model 9
m9<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen) +
          lag(local) +
          lag(state)*lag(goals_ratio) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi1,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m9)

#####Output####
texreg(list(m5,m6,m7,m8,m9), 
       digits = 3,
       custom.model.names = c("(5)","(6)","(7)","(8)","(9)"),
       custom.coef.map = list("lag(cgr)" = "CGR",
                              "lag(citizen)" = "Citizen group",
                              "lag(local)" = "Local gov",
                              "lag(state)" = "State gov",
                              "lag(goals_ratio)" = "Goals",
                              "lag(cgr):lag(goals_ratio)" = "CGR x goals",
                              "lag(citizen):lag(goals_ratio)" = "Citizen group x goals",
                              "lag(local):lag(goals_ratio)" = "Local gov x goals",
                              "lag(state):lag(goals_ratio)" = "State gov x goals",
                              "lag(representation)" = "Representation",
                              "lag(representation^2)" = "Representation sqr",
                              "lag(project)" = "N-Project",
                              "lag(project^2)" = "N-Project sqr",
                              "lag(kind_ratio)" = "In-kind share",
                              "forest" = "Forest land use (%)",
                              "air.temp" = "Temperature (°C)",
                              "ppt" = "Precipitation (mm)",
                              "w_wqi1" = "Spatial lag",
                              "(Intercept)" = "Constant"))

####Table A3####
vars <- names(df)[-(c(1:4, (ncol(df)-1):ncol(df)))]
des_stat<- df %>%
  dplyr::select(vars) %>%
  map_df(~get_summary_stats(.x), .id = "variable")

des_table<-xtable(des_stat, digits=2)
print.xtable(des_table)

####Table A4####
moran_df <- data.frame(
  year = numeric(),
  moran_I = numeric(),
  expected_I = numeric(),
  z_score = numeric(),
  p_value = numeric()
)

years <- unique(df$year)

for(year in years) {
  df_year <- df[df$year == year,]
  df_year <- na.omit(df_year)
  
  if(nrow(df_year) == 0) {
    print(paste("No observations in year", year))
    next
  }
  
  df_year_coords <- df_year[,c("lon","lat")]
  
  dist_matrix <- sp::spDists(as.matrix(df_year_coords))
  inv_dist_matrix <- 1 / dist_matrix
  diag(inv_dist_matrix) <- 0
  
  listw <- mat2listw(inv_dist_matrix, style="W")
  
  test <- try(moran.test(df_year$wqi, listw), silent = TRUE)
  
  if(class(test) == "try-error") {
    print(paste("Error in year", year))
    next
  }
  
  moran_I <- as.numeric(test$estimate["Moran I statistic"])
  expected_I <- as.numeric(test$estimate["Expectation"])
  z_score <- as.numeric(test$statistic)
  p_value <- as.numeric(test$p.value)
  
  moran_df <- rbind(moran_df, data.frame(year, moran_I, expected_I, z_score, p_value))
}
moran_df <- moran_df[order(moran_df$year), ]

# table
print(moran_df)
moran_table<-xtable(moran_df, digits=2)
print.xtable(moran_table)

####Table A6####
# Longitudinal Analysis with the Distance-based Spatial Lag
# Model 10
m10<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen) +
          lag(local) +
          lag(state) +
          lag(goals_ratio) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi2,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m10)

# Model 11
m11<-plm(wqi ~ 
          lag(cgr)*lag(goals_ratio) +
          lag(citizen) +
          lag(local) +
          lag(state) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi2,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m11)

# Model 12
m12<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen)*lag(goals_ratio) +
          lag(local) +
          lag(state) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi2,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m12)

# Model 13
m13<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen) +
          lag(local)*lag(goals_ratio) +
          lag(state) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi2,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m13)

# Model 14
m14<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen) +
          lag(local) +
          lag(state)*lag(goals_ratio) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi2,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m14)

#####Output####
texreg(list(m10,m11,m12,m13,m14), 
       digits = 3,
       custom.model.names = c("(10)","(11)","(12)","(13)","(14)"),
       custom.coef.map = list("lag(cgr)" = "CGR",
                              "lag(citizen)" = "Citizen group",
                              "lag(local)" = "Local gov",
                              "lag(state)" = "State gov",
                              "lag(goals_ratio)" = "Goals",
                              "lag(cgr):lag(goals_ratio)" = "CGR x goals",
                              "lag(citizen):lag(goals_ratio)" = "Citizen group x goals",
                              "lag(local):lag(goals_ratio)" = "Local gov x goals",
                              "lag(state):lag(goals_ratio)" = "State gov x goals",
                              "lag(representation)" = "Representation",
                              "lag(representation^2)" = "Representation sqr",
                              "lag(project)" = "N-Project",
                              "lag(project^2)" = "N-Project sqr",
                              "lag(kind_ratio)" = "In-kind share",
                              "forest" = "Forest land use (%)",
                              "air.temp" = "Temperature (°C)",
                              "ppt" = "Precipitation (mm)",
                              "w_wqi2" = "Spatial lag",
                              "(Intercept)" = "Constant"))

####Table A8####
# Model 15
m15<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen) +
          lag(local) +
          lag(state) +
          lag(goals_ratio) +
          lag(align_ratio_per) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi1,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m15)

# Model 16
m16<-plm(wqi ~ 
          lag(cgr)*lag(goals_ratio) +
          lag(citizen) +
          lag(local) +
          lag(state) +
          lag(align_ratio_per) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi1,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m16)

# Model 17
m17<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen)*lag(goals_ratio) +
          lag(local) +
          lag(state) +
          lag(align_ratio_per) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi1,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m17)

# Model 8
m18<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen) +
          lag(local)*lag(goals_ratio) +
          lag(state) +
          lag(align_ratio_per) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi1,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m18)

# Model 19
m19<-plm(wqi ~ 
          lag(cgr) +
          lag(citizen) +
          lag(local) +
          lag(state)*lag(goals_ratio) +
          lag(align_ratio_per) +
          lag(representation) +
          lag(representation^2) +
          lag(project) +
          lag(project^2) +
          lag(kind_ratio) +
          forest +
          air.temp +
          ppt +
          w_wqi1,
        data = df,
        index = c("huc8", "year"),
        effect = "twoways",
        model = "random")
summary(m19)

#####Output####
texreg(list(m15,m16,m17,m18,m19), 
       digits = 3,
       custom.model.names = c("(15)","(16)","(17)","(18)","(19)"),
       custom.coef.map = list("lag(cgr)" = "CGR",
                              "lag(citizen)" = "Citizen group",
                              "lag(local)" = "Local gov",
                              "lag(state)" = "State gov",
                              "lag(goals_ratio)" = "Goals",
                              "lag(cgr):lag(goals_ratio)" = "CGR x goals",
                              "lag(citizen):lag(goals_ratio)" = "Citizen group x goals",
                              "lag(local):lag(goals_ratio)" = "Local gov x goals",
                              "lag(state):lag(goals_ratio)" = "State gov x goals",
                              "lag(align_ratio_per)" = "Goal alignment",
                              "lag(representation)" = "Representation",
                              "lag(representation^2)" = "Representation sqr",
                              "lag(project)" = "N-Project",
                              "lag(project^2)" = "N-Project sqr",
                              "lag(kind_ratio)" = "In-kind share",
                              "forest" = "Forest land use (%)",
                              "air.temp" = "Temperature (°C)",
                              "ppt" = "Precipitation (mm)",
                              "w_wqi1" = "Spatial lag",
                              "(Intercept)" = "Constant"))

